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Abstract 

Given two metastable states A and B of a biomolecular system, the problem is to 
calculate the likely paths of the transition from A to B. Such a calculation is more 
informative and more manageable if done for a reduced set of collective variables chosen 
so that paths cluster in collective variable space. The computational task becomes 
that of computing the "center" of such a cluster. A good way to define the center 
employs the concept of a committor, whose value at a point in collective variable 
space is the probability that a trajectory at that point will reach B before A. The 
committor "foliates" the transition region into a set of isocommittors. The maximum 
flux transition path is defined as a path that crosses each isocommittor at a point 
which (locally) has the highest crossing rate of distinct reactive trajectories. (This 
path is different from that of the MaxFlux method of Huo and Straub.) It is argued 
that such a path is nearer to an ideal path than others that have been proposed with 
the possible exception of the finite-temperature string method path. To make the 
calculation tractable, three approximations are introduced, yielding a path that is the 
solution of a nonsingular two-point boundary-value problem. For such a problem, one 
can construct a simple and robust algorithm. One such algorithm and its performance 
is discussed. 



1 Summary 

Considered here is the problem of computing transition paths of conformational change, 
given two different metastable states of a biomolecule. One motivation for this is to facilitate 
the accurate calculation of free energy differences. Another motivation is to determine the 
existence and structure of transition states and intermediate metastable states. The latter are 
possible targets for inhibitors of enhanced specificity in cases where a family of proteins have 
active sites with very similar structure. A good example of this situation is the Src tyrosine 
kinase family [39], which has long been implicated in the development of cancer. For this 
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system there are already computational results [TTJ |2H EE] , supported by experiment [25J, 
for the transition path from an active catalytic domain to an inactive catalytic domain. 

Some approaches to this problem generate ensembles of trajectories based on the equa- 
tions of motion. Notable examples are transition path sampling [2] and Markov state mod- 
els [32]. Applying such methods to large proteins (without compromise) would appear to re- 
quire exceptional computing capabilities, so here we pursue a more theoretical approach that 
avoids "direct numerical simulation." Such approaches, like those in Refs. [131 E21 EZ] , seek 
to characterize one (or several isolated) "representative" reaction paths connecting two given 
metastable states, each path representing a bundle or cluster of trajectories. Here we adopt 
a well developed and tested theory, namely, transition path theory (TPT) [8j [201 ED E3] . 
Additional references on computing transition paths are found in Ref. [26J. In general, it 
may also be of interest to calculate (i) the reaction rate for each bundle, or, at least, the 
relative rate for different bundles, and (ii) the potential of mean force. Here we consider only 
the calculation of the path itself. 

In a nutshell, this article embraces a certain aspect of TPT and carries it to a logical 
conclusion, obtaining a formula, an implementation, and a proof of concept. The claim 
is that we can compute a path that is closer to the ideal than the minimum free energy 
path (MFEP) [20] and that, in a couple of respects, is better than the path of the finite- 
temperature string (FTS) method [291 EE], though inferior in another (important) respect. 
Additionally, the formula for the path is computationally more attractive than the formula 
that underlies either the path of the FTS method or the MFEP. 

1.1 Outline and discussion 

There are two distinct steps in getting a solution: The first is to define the problem without 
concern for the methods to be employed (other than taking into account the intrinsic difficulty 
of the problem). Defining a problem apart from a method gives a more concise definition. 
Also, by not guessing about what is feasible computationally, one may avoid unnecessary 
compromises. The second step is to construct a method and algorithm. 

Given two metastable states A and B of a biomolecular system, the aim is to calculate 
the likely paths of the transition from A to B. Such a calculation is more informative and 
more manageable if done for a reduced set of collective variables, functions of the system 
configuration x, 

Ci = £iO), C2 = &0*0 > • • • , Cu = £v(x), abbreviated as ( = £(x), 

chosen so that paths cluster in collective variable space. The computational task becomes 
that of computing the "center" of such a cluster. A good way to define the center employs the 
concept of a committor, whose value at a point in collective variable space is the probability 
that a trajectory at that point will reach B before A. The committor "foliates" the transition 
region into a set of committor isosurfaces known as isocommittors. The maximum flux 
transition path (MFTP) is defined as a path that intersects each isocommittor at a point 
which (locally) has the highest crossing rate of distinct reactive trajectories. The MFTP is 
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not to be confused with MaxFlux method [HUH]; it differs in several respects, in particular, 
the MFTP considers the flux of only those trajectories that are reactive (by using a result 
from TPT). A more detailed account of the problem definition is given in Section [2j 

The minimum free energy path has been used for some time to represent reactive trajec- 
tories in collective variable space. Only fairly recently has its relationship to reactive trajec- 
tories been explained. The article [20] applies large deviation theory to show that the MFEP 
is the most probable path in the zero temperature limit of dynamics on a free energy surface 
defined at finite temperature. Hence, the MFEP is an inherently inconsistent construct and 
it is useful only to the extent that it represents fully finite-temperature trajectories. In fact, 
it does this fairly well on the simple tests reported here. Other fully finite-temperature con- 
structs have been proposed: the finite-temperature string method in collective variable space 
(Sec. IV. B of Ref. [2H]) and the swarm-of-trajectories string method [25], which constructs 
a Brownian dynamics model on the fly and constructs a path whose tangent is the most 
probable direction. How they differ from the MFTP is detailed in Section [2} 

To make the calculation tractable, three approximations are introduced. To make the 
committor a more accessible quantity, the set of paths is approximated by a Brownian dy- 
namics model, resulting in a boundary value problem in //-dimensional space. Then the 
number of space dimensions is reduced to one by assuming most of the transition paths are 
contained in a tube, resulting in a two-point boundary-value problem with 2v unknowns. A 
third approximation reduces this to v unknowns, whose solution is a maximum flux tran- 
sition path. The resulting equations involve a free energy gradient term and an explic- 
itly temperature-dependent curvature term. Specifically, the maximum flux transition path 
( = Z(s), < s < 1, is defined by the condition that 



holds for ( = Z(s) where (3 is the inverse temperature, F(() is the free energy profile, 
D(() is a proto-diffusion tensor depending on masses and £, and the subscript s denotes 
differentiation (d/ds). In the high temperature limit, the path becomes a straight line. 
In the low temperature limit, the path becomes an MFEP. At zero temperature the path 
will have cusps at some intermediate local minima, which presents difficulties if free energy 
profiles or relative reaction rates are to be determined. This formula is a key result of this 
article. Details are given in Section [3j 

The temperature-dependent curvature term not only provides a finite temperature cor- 
rection to the MFEP, but it yields a nonsingular second order ordinary differential equation, 
amenable to standard techniques — except for the need to do computationally intensive sam- 
pling to evaluate terms in the differential equation. An existing set of algorithms for the 
MFEP [61 [20] applies equally well to the MFTP. The equation is discretized using upwinded 
differencing and solved using the semi-implicit simplified string method [33]. (A notable 
alternative is the nudged elastic band method, introduced in Ref. [2].) Algorithmic details 
are provided in Section |4j 

Section [5] compares the MFTP to the MFEP on numerical examples. First, an artificial 
problem in full configuration space is solved to demonstrate the effect of the curvature term 
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of the MFTP. (A problem in full configuration space is equivalent to a problem in collective 
variable space with perfect sampling.) In particular, the necessity of using an adaptive mesh 
for the MFEP is demonstrated. Then alanine dipeptide in vacuum is solved using the <ft, ip 
dihedrals as collective variables. For the transition path from CV ax to CV eq as in Ref. [20J, 
the computational cost for calculating the MFTP and the MFEP is almost same. However, 
for a transition path from CV cq to C' 7cq through CV ax shown in Ref. [2H], the MFEP has a 
cusp at C 7ax and the computational cost for finding such a cusp is expensive. On the other 
hand, the MFTP smooths out the cusp and the computational cost is reduced. 

An open source implementation of the MFTP method is available [ID] as a relatively sim- 
ple set of Python modules with examples using pure Python, CHARMM [I] , and NAMD [28J . 

1.2 Conclusions 

For alanine dipeptide, the MFEP, MFTP, and FTS method paths are quite similar. On a 
contrived problem with a rough energy landscape, e.g., Figure 2 in Ref. [36], the FTS method 
path gives a much better result. On a different contrived problem given in Section [5TT| the 
MFTP gives a much better result. Contrived examples are relevant because computational 
techniques are sometimes applied in extreme situations for which they may not have been 
designed. In terms of quality, the MFTP ranks higher than the MFEP but lower than the 
FTS method path (because the latter addresses the more serious difficulty of multiple local 
minima) . 

The minimum free energy path (and that of the FTS method) can have cusps at some 
intermediate metastable states, which makes it unsuitable for defining an isocommittor, un- 
suitable for defining a reaction coordinate, and harder to compute. Computational difficulties 
include the need for an adaptive mesh and a greater number of iterations until convergence. 

2 What is the problem? 

We begin by defining an ensemble of transition paths from A to B: For simplicity, assume 
the molecular system obeys Newtonian dynamics with potential energy function U (x) and a 
diagonal matrix M of atomic masses. Positions x and momenta p satisfy x = X(t), p = P(t) 
where (d/dt)X(t) = M~ l P{t) and (d/dt)P(t) = -VU(X{t)). Initial values are drawn from 
a Boltzmann-Gibbs distribution p(x,p): positions x from probability density const e~ T ^ 
and momenta p from a Maxwell distribution. Imagine an extremely long trajectory. The 
trajectory enters and leaves A and B many times yielding a huge set of reactive paths from 
A to B. (A reactive path is a piece of the trajectory outside of A and B that comes from A 
and goes to B.) 

Generating an ensemble of trajectories is extremely demanding computationally. And, 
even if this were possible, what would the user do with all the data? By answering such 
a question, we might well avoid the task of computing trajectories. It is likely that one 
would cluster the trajectories to produce a concise description. Therefore, one might instead 
directly determine such a concise description. Specifically, if the paths cluster into one or 
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several distinct isolated bundles/tubes/channels/pathways, one might compute a "represen- 
tative path" for each cluster. This idea is developed in the paragraphs that follow. 

However, transition paths might not cluster adequately — in full configuration space. As- 
sume, though, there is a smaller set of collective variables, ( = such that in £-space, 
paths cluster into one or several distinct isolated channels connecting two separated subsets 

and of collective variable space. Otherwise, there is little of interest to compute. A 
typical example of collective variables is <f)/ip angles along a peptide backbone. Once the 
collective variables are specified, the problem is to calculate a path in collective variable 
space, ( = Z(s), < s < 1, connecting to where the transition paths are concen- 
trated. Along with a parameterization of the path in collective variable coordinates, would 
be a realization of it in cartesian coordinates, so once the path is generated, structures can 
be studied as well. A drawback of this approach is the need to identify an appropriate set 
of collective variables. Indeed, defining suitable collective variables is an important research 
problem [T7]. 

We want a minimal set of collective variables subject to two conditions: First, the coor- 
dinates C must suffice to describe states A^, B^ in £-space corresponding to A, B. Second, 
coordinates ( must also be rich enough to "express the mechanism of conformational change" 
along the transition path. To make the second condition more precise, we introduce the no- 
tion of "quasi-committor." 

To measure the progress of a transition, there is a natural reaction coordinate, known as 
the committor. This concept of a commitment probability was introduced by Onsager |23j, 
and the abbreviated term "committor" was introduced in Ref. [3] (p. 9236), which they 
defined as follows: For each point x in configuration space, consider a trajectory starting 
with X(0) = x and velocities drawn at random from a Maxwell distribution, and define the 
committor q(x) to be the probability of reaching B before A. Since it is the coordinates of 
the collective variables that are of interest, it is natural also to define a quasi-committor: 
For each point (, consider a trajectory starting with random initial values conditioned on 
£(x) = ( and define the quasi-committor q(() to be the probability of reaching B^ before A%: 

q(Q = Pr(£(X(t)) reaches B^ before | £(X(0)) = C)- 

We could say that the variables ( = are rich enough to express the mechanism 
of conformational change if the quasi-committor q(() has no local minima or maxima out- 
side of A^ and B^ (except for regions of negligible probability). Otherwise, there is some 
unexpressed degree of freedom important to the transition. As an example, suppose that 
virtually all trajectories stay within a narrow tube having a geometry in full configuration 
space illustrated by [TJ Suppose that the free energy profile as a function of arc length along 
the transition tube is much higher in the backward section than it is in the two forward sec- 
tions. Then most of the increase in the quasi-committor as a function of arc length occurs 
in the middle section of the tube. Consequently, the variation in the quasi-committor, as a 
function of the ill-chosen collective variable ( corresponding to the horizontal axis, will be 
dominated by this middle section of the tube. This results in a graph of q(() that increases 
at the beginning and end of its range but decreases in the middle part. In addition to q(() 
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having no local extrema, it is desirable that q(£(x)) ~ q(x). The quality of the collective 
variables can be checked in principle by calculating quasi-committor values q(Q at points 
along the path from dynamics trajectories. 



Figure 1: A schematic illustration of a poor choice of collective variables. The horizontal 
axis is collective variables, and the vertical axis is unrepresented degrees of freedom. The 
collective variables fail to indicate the progress of the reaction. 

Two approaches have been proposed for defining the center of a cluster of paths in £-space: 

(i) a most probable path, e.g., a swarm-of-trajectories string method [26], El] path, and 

(ii) a path that intersects each isosurface of the quasi-committor at a center of the collection 
of points where reactive trajectories cross that isosurface, e.g., a finite temperature 
string method [29] path and a maximum flux transition path. 

An MFEP is a limiting case of both approaches. (The MFEP is obtained from these various 
approaches by letting j3 — > oo in the path formula but not in the definition of the free energy 
profile.) Defining the objective is a compromise between (i) best capturing the object of 
interest and (ii) simplicity. 

One problem with seeking the most probable path is that it is unclear how to assign 
relative probabilities to paths. More importantly, the most probable path tends to be a path 
of minimum energy, and it is not clear — a priori — that this is a "representative" path. For 
Hamiltonian dynamics, it would seem that the probability that we attach to a path would 
be proportional to exp(— (3E) where E is the energy. Hence, the most probable path is 
the one with just enough energy to surmount the potential energy barriers. For stochastic 
dynamics, the explanation of how to assign probability to paths is quite complicated — if 
paths of different durations are being compared. An explanation for Brownian dynamics is 
possible using Freidlin-Wentzell theory and the assumption of vanishingly small noise (see 
Appendix A of Ref. |20j). It is reassuring though that the results of Freidlin-Wentzell theory 
agree with those of TPT in the zero-temperature limit (for F(Q held fixed). 

For defining a path in terms of an intersecting point on each isosurface of a quasi- 
committor, one needs 

(i) a definition for the distribution of crossing points of reactive trajectories through a 
quasi-commitor isosurface and 
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(ii) a definition of centrality, e.g., mode, median, or mean. 

We consider each of these in turn. 

The finite-temperature string method defines the distribution of crossing points of re- 
active trajectories in a way that includes recrossings. A subsequent article illustrates 
the dramatic distortions that arise by including recrossings, and it emphasizes crossings of 
a surface by distinct reactive trajectories instead of all crossings by reactive trajectories. 
They define such a distribution in terms of the net crossings of reactive trajectories across 
each infinitesimal piece of a surface. It is not obvious, however, that this necessarily gives 
nonnegative values on the isosurface of a quasi-committor, so, instead, we use the density of 
last crossings by reactive trajectories, called last hitting points in Ref. [8]. For the Brownian 
dynamics approximation developed in the next section, these two measures are identical. 

Consider now the question of defining the center. Let j(Q denote the density associated 
with a definition for the distribution of crossing points of reactive trajectories through a quasi- 
committor isosurface. One choice for the center is the point of highest probability, In other 
words, seek the path ( = Z(s), < s < 1, each of whose points Z(s) is a local maximum of 
the density on the quasi-committor isosurface £ passing through Z(s). This is what we 
use for the MFTP. Another choice, associated with the finite-temperature string method, is to 
construct the path from the mean value (' on each quasi-committor isosurface E: the point 
(' that minimizes J s |£' — C| 2 j(C)dC Although this notion is a superior measure of centrality, 
it is more complicated to explain. In practice, methods for finding a maximum are designed 
only to find a local maximum, which is what we do for the MFTP. This is satisfactory if there 
is a choice of collective variables that produces a free energy landscape free of roughness at 
the scale of the thermal energy [36]. In any case, the equations defining a center-of-density 
path are intrinsically more expensive computationally to solve than those for the MFTP, 
because they require averaging on quasi-committor isosurfaces q(Q = constant (in addition 
to conditional averages on collective variable isosurfaces £(:r) = £ in full configuration space) 
rather than merely determining a (local) maximum. 



3 A method 

As stated previously, computing q(() is not feasible. Consequently, we derive a method, 
which employs three uncontrolled approximations — a controlled approximation being one 



that can be made arbitrarily accurate with sufficient computational effort. Subsection 3.1 



approximates paths in collective variable space by those of Brownian dynamics; Subsec- 



tion 3.5 assumes most paths lie in a tube where isocommittors are planar; and Subsection 3.6 
assumes that on average the trajectories are parallel to the path. The basic ingredients of 
much of this development are present in the literature but scattered among several articles. 
Here they are combined to produce equations from which we derive the MFTP. 
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3.1 Brownian dynamics approximation of collective variable paths 

The probability density function (p.d.f.) for £(x) is 



Pd0 = " 0> = J J ~ C)p(x, P )dxdp 

where S(Q = S((x)5((2) • ■ ■ S(Cv)- Let (-)^ be the expectation for the conditional density 
p(x,p\£(x) = C): 

(<Kf W - 0) 

In Appendix [A] is an adaptation of an argument from Ref. [20] (Sec. Ill, A and B) 
suggesting that as an approximation to q(C), we should seek a function q(Q that minimizes 
a certain functional I(q) that can be expressed in terms of collective variables £. Define the 
free energy F(() for coordinates £ = by 

con^e-^=p € (C) = (^)-C)). (1) 
Also define a proto-diffusion tensor D by 

0(0 = ir 1 <6 B (*)M- i e e (x) T > e . 

(There is freedom in the scaling of D. We use this freedom to make Eq. Q agree with 
an alternative derivation of the Brownian dynamics, in which one assumes instantaneous 
relaxation of the degrees of freedom not represented by the collective variables. The tensor 
-D(C) fails to be a diffusion tensor because it is missing a time scale factor.) The functional 
is then 

I(q) = consts j e-^Vg(C) T J D(C)Vg(C)dC (2) 

where the integral is over the transition region outside of and subject to q(() = on 
the boundary of and q(() = 1 on the boundary of B^. 

The corresponding Euler-Lagrange equation for q(Q is the Smoluchowski (backward Kol- 
mogorov) equation: 

- V • e-P F ®D(Q Vg(C) = 0, (3) 

subject to q(() = on the boundary of A^ and q(() = 1 on the boundary of B^. 

The function q that satisfies the Smoluchowski equation subject to the given boundary 
conditions can be shown to be the exact committor function for paths C — C( r ) m collective 
variable space generated by the Brownian dynamics 

= -WC)VF(C) + (V ■ D(C)) T + V2D 1/2 (() V (t) (4) 

where Di/ 2 Dj^ 2 = D and r/(r) is a collection of standard white noise processes. The fact that 
t is an artificial time does not affect the committor. In principle, the assumption q(() q(() 
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can be checked a posteriori by comparing committor values of the Brownian dynamics to 
the quasi-committor values of actual dynamics. 

Reference [20] (Sec. III.C) appears to suggest that the Smoluchowski equation uniquely 
specifies dynamics except for scaling of time: If the Smoluchowski equation ^ is satisfied by 
committors q(Q for arbitrary sets A'^ and B'^ in collective variable space, then trajectories 
whose committor functions satisfy Eq. ^ must have paths that are those of the Brown- 
ian dynamics. Hence, paths in collective variable space can be generated with the proper 
probabilities from the system of stochastic differential equations. 

3.2 Last hitting-point distribution 

Appendix [B] considers the rate at which reactive trajectories cross an arbitrary surface £ that 
separates collective variable space into two parts, one containing and the other containing 
B^. The result given there is that the rate of the last crossing of E by reactive trajectories 
is given by the integral 

/ J(C) • h(QdS 
where n{C) points to the side containing B^, and 

AO = P 5 (C)£(C)Vg(C) 

is the last hitting-point flux. The choice of the last hitting point to represent the point 
where a reactive trajectory crosses an isocommittor is somewhat arbitrary. Therefore, it is 
gratifying to know that the expression for J(Q also gives the net flux and the first hitting- 
point flux of reactive trajectories. 

The normal to an isocommittor is given by h(() = Vq(Q/\Vq(Q\, so the distribution of 
last hitting points on an isocommittor is proportional to 

3(0 = ^(C)Vg(C) T ^(C)Vg(C)/|Vg(C)|. 

3.3 Defining the path 

For computation it is convenient to label the isocommittors with the path parameter s. In 
particular, denote by S(s), < s < 1, the isocommittor passing through ( = Z(s). Write 
<?( s ) = <?(-^( s )) an d define <r(() implicitly by 

q(0 = <?H0). (5) 

In this way the committor q(() is decomposed into two independent parts: one part a(Q 
specifies the isocommittor label and the other part q(s) calibrates the isocommittors. On an 
isocommittor the gradient Vg(C) = ?s(c"(C))Vcr(C), so the normal flux is 

J(0 = ^(^(C))p ? (C)Va(C) T J D(C)Va(C)/|Va(C)| (6) 
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(recalling that the subscript s denotes differentiation d/ds). Note that q s contributes only a 
scale factor to j(C), so the center of intensity on E(s) does not depend on q s . 

Each point Z(s) on the desired path maximizes the last hitting-point flux j(Q on the 
isocommittor q(Q = q(Z(s)). Hence, Vj(Z(s)) || Vq(Z(s)). To keep the derivation indepen- 
dent of the calibration q(s), introduce a vector n(s), not necessarily normalized, such that 
n(s) || Vq(Z(s)). Hence, 



3.4 The localized tube assumption 

Assume there exists a tube connecting to such that (i) on each isocommittor, regions 
of high j{C) are concentrated in the tube, (ii) each isocommittor is nearly planar in the tube, 
and (iii) -D(C) is nearly constant on each isocommittor within the tube. This scenario is 
illustrated in [2] below. 



Figure 2: Shading indicates contours of free energy, thin curves denote isocommittors, 
ellipses enclose concentrations of crossing points from reactive trajectories, and the thick 
curve is the center. 

Exploit the localized tube assumption by approximating the isocommittor through Z(s) 
as a plane n(s) with normal n(s). Hence, the isocommittor surface S(s) : a(() = s has the 
simple description of a hyperplane, 




(7) 







10 



Also, approximate D(() by D(a(()) where D(s) = f D(Z(s)). These approximations (see 
Ref. [33] (Sec. 6.6.1)) are sufficient to define a practical method (see Ref. [7] (Sec. 12)). 
The unknown direction vector n(s) is to be chosen to minimize the integral I(q) of Eq. ^ 
restricted to some tube. For simplicity the boundary points Z(0) and Z(l) can be moved 
to points in and that locally minimize F((). In this way the problem of solving for 
a committor of many variables is reduced to that of a one-dimensional calculation along the 
length of the tube. 

It remains to derive the condition that determines Z(s). This is done in Appendix [Cj 
where it is shown that the condition is 

-^> + ^""- 

3.5 The maximum flux transition path 

Although the localized tube assumption is sufficient for defining a practical method, the 
method would not be simple, so we make an additional simplifying assumption: Assume the 
flux J(() points in the direction of the path so that J(Z(s)) \\ Z s (s) or D(Z(s))Vq(Z(s)) \\ Z s (s) 
whence 

n(s) || D(Z(s))' 1 Z s (s). 
The result is a maximum flux transition path 

- F?F{Z) + { ^(Z) Z 'z a 11 °W lZ °- ^ 

(The simplifying assumption is justified, for example, if the probability is strongly peaked 
around the path, resulting in most of the probability contained in a narrow tube with a flux 
J(C) pointing in the direction of the tube and the path.) This assumption is also made for 
the FTS method, see Eq. (14) of Ref. [23] and Sec. II. A. of Ref. [3S]. Geometrically, this 
condition means that instead of having the free energy gradient vanish orthogonal to the 
path, it is balanced by a "centripetal" force, which reduces curvature and avoids cusps. 

To express condition ^ as an equation, write it as -P(ZjD~ 1 Z s )DVF + D{p- x Z s ) s = 
XZ S where A is a scalar and premultiply by Zj to obtain an expression for A. After eliminating 
A, the equation becomes 

(/ - R)D(D~ 1 Z S ) S = (I - \l)P{Z T s D- 1 Z s )DVF (10) 

where the projector 

U = Z s Z]/(ZjZ s ). 



Note that, if D is constant, the limit (3 — > for Eq. (10) gives a geodesic Z ss = 0, which 
is the desired result. 

In the two-dimensional case with D = I, the Euclidean length of (I— U)D(D^ 1 Z S ) S / (ZjD~ 1 Z s 
is exactly equal to the curvature, which is defined to be the reciprocal of the radius of cur- 
vature. To see this, note that this is true if we parameterize with (actual) arc length and 
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note also that the curvature term is independent of parameterization (which can be checked 
analytically). 

If we normalize the parameterization using (ZjZ s ) s = 0, this implies ZjZ ss = and 
IiZ ss = 0. Adding this last equation to Eq. ( [To] ) gives a nonsingular second order ordinary 
differential equation for Z(s): 

Z ss = (I - II) ^{Z T S D- 1 Z S )DVF + D S D- X Z S ) . 

Values obtained from constructing the path can be used to calculate the free energy 
F(Z(s)) along the path, 

F(Z(s)) - F(Z(0)) = f\F{Z{s')YZ s {s')ds'. 

Jo 

However, F(Z(s)) is not a potential of mean force for the transition. 
3.6 The minimum free energy path 

The simplifying assumption of the preceding subsection, which is used to derive the MFTP, 
is valid in the limit /3 — > oo in the Brownian dynamics approximation; see Ref. [33] (Sec. 6.6) 
and Ref. [20] (App. A). A more systematic derivation might therefore neglect the curvature 
term. The result would be a minimum free energy path 

Z s || -(3D(Z)VF(Z). 

Each point ( = Z(s) on the MFEP is a local minimum of F(() in the hyper-plane orthogonal 
to D(Z(s))~ l Z s {s). 

One difference from an MFTP is that an MFEP can have a cusp at an intermediate 
local minimum. If the path passes sufficiently close to a local minimum ( = (o of F(Q, 
then for a short section of the path, ( = Z(s), a < s < b, a quadratic approximation to 
F(() is accurate. Assume D = constant and F(() = |(£ — ( ) T A(( — Co)+ constant, where 
A is symmetric positive definite. The MFEP is then defined by Z s \\ — (3DA(Z — £ ). 
Perform a change of variables, Y = /3~ 1 ^ 2 Q T D^ 2 (Z — £o) where QAQ T is a diagonalization 
of D~[j 2 ADi/2- The MFEP for Y(s) is hence given by Y s || — AY. For simplicity, suppose 
that Y = [x,y] T , that x(a) < < x(b), and that A = diag(A,/i) with A > /i. The path is 
hence defined by y s /(^y) = x s /(Xx), which can be integrated to yield the path 

_ f (x/x(a)Y /x y(a), x(a) <x<0, 
y ~\ (x/x(b)y/ x y(b), 0<x< x(b), 

which has a cusp at x = 0. 

The FTS method path is also likely to suffer from the presence of cusps, because for a 
harmonic potential, the average position is the same as the most probable position. 

The presence of cusps undermines the localized tube assumption. In particular, the 
assumption of isocommittors being approximately planar breaks down at a cusp. This poses 
a difficulty when computing transition rates, which depends on existence of isocommittors. 
Additionally, cusps complicate the numerical approximation of paths. 
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4 An algorithm 



An algorithm for calculating a transition path employs a progression of four controlled ap- 
proximations: discretization of the path £ = Z(s) and the equations that define it; a finite 
number of iterations for the solution of nonlinear discrete equations; use of restraints for 
constrained sampling; and finite sampling. 

4.1 Discretization 

The path Z(s), < s < 1, is approximated as a piecewise polynomial with break points 
= s < s i < • • • < s j — 1- Here we choose a uniform mesh s = 0, As, . . . , 1 and obtain 
the path by piecewise linear interpolation. Thus the problem is reduced to determining 
unknown nodal values Zj ~ Z(sj), j = 0, 1, . . . , J, each representing a replica of the system 
in a different configuration. 

It is convenient for computation to use for the path parameter s the arc length along the 
path divided by the total length of the path. In such a case, |Z s (s)| is constant. The arc 
length normalization becomes 

\Z j+ i - Zj\/As = \Zj - Z^l/As, j = 1, 2, . . . , J - 1. 

Condition ^ is written as 

-PDVF + -Z ss - -D a D~ x Z s || Z s 

c c 

where c = ZjD~ 1 Z s . 

This is discretized by the finite difference scheme 

(Zs), || 9j , where 9j d = f -^(VF), + i Z ' +1 - + ^ - -^ZT 1 ^ 

and where 

Cj = ^As- 2 (A_Z7 J D7 1 A_Z J + A + Z7 J D7 1 A + Z J ), (11) 
(DsD^Zs), = ^As- 2 (A_ J D jJ D7 1 A_Z J + A +J D iJ D7 1 A+Z J ), (12) 

with 

A±L>j = T(Dj - D i±1 ), and A ± ^ = t(Z 3 - Z j±1 ). 
We choose upwinded differencing for (Z s )j based on the direction of the modified mean force 

gr- 




in the unlikely event that both conditions are satisfied, the choice is dictated by the arc 
length normalization step of the simplified string method to be discussed next. 

For the MFEP, cusps can occur at some intermediate local minima, requiring an adaptive 
mesh to resolve. 
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4.2 Solution of nonlinear discrete equations 



A second component of the algorithm is an iterative method for achieving rapid local con- 
vergence given a plausible initial guess. 

Because of its simplicity and demonstrated effectiveness, we adopt the semi- implicit sim- 
plified string method used in Ref. [35] (Eq. (11)). To determine a path, begin with an initial 
guess and generate successive improvements by alternating between moving the points of the 
curve Zj in the direction gj given by condition ^ and reparameterizing. 

The first step of each iteration is to solve the following equations for the Z*: 



Z*-Z, 



1 Z* - 2Z* + Z* 1 

- 3+1 a \ 3 -=±- -{D.D- 1 Z.) j -PD j &F i ) j , j = l,2,-. ,J-1, 



T 2 Cj As 2 Cj 



T 2 



-PDjiVF)^ j = 0,J, 



where Cj and (D s D~ 1 Z s )j are given in Eqs. (11) and (12). (The extra factor r provides the 
time scale factor missing from D.) 

Then the normalization adjustment is to choose the {Zj} to be equidistant along the 
resulting curve: 

si = 0,s* = s*_ 1 + \Z*-Zl 1 \, 
Z*(s) = piecewise linear interpolation of {(s*/s*j,Z*)}, < s < 1, 

Z new = Z*(j/J). 

It can be shown that if the semi-implicit simplified string method converges, the resulting 
points Zj satisfy a nonstandard discretization of the differential equation containing r as a 
parameter. In the limit r — > 0, the discretization becomes upwinded differencing. 

For large systems, targeted molecular dynamics [31] has been used to get an initial 
path [11] [12] . Another potentially promising but quite different approach is rigidity analy- 
sis 



4.3 Conditional averages 

Evaluation of VF and D at break points involves sampling on hyper-surfaces {x : £(a:) = Zj} 
of configuration space. 

For calculating such conditional expectations, the Dirac delta function S(s) can be ap- 
proximated by the p.d.f. of a Gaussian 5 £ (s) = {2ne 2 )^ 1 ^ 2 exp(— s 2 /(2e 2 )). Note 

S e (£( x ) - C)e-^ = (27re 2 r»/ 2 e- mx ^ 

where 

v 1 

U(x;C) = U(x) + ^ Ui (x,(i), and Ui (x,Ci) = ^glMx) - Ci) 2 - (14) 
i=i P 
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Then, (0(x))^ = (0(x)6 e (£(x) — ()) / (5 £ (t;(x) ~~ 0) is nothing but an average using U(x;Q. 
The effect is that of using restraining potentials instead of constraints. These restraints 
should be as strong as possible without restricting the step size used in the sampling. From 
const^ exp(— f3F(()) = (5 E (£(x) — £)}, we have 

vF(c) = -^<e(*)-0c- 

4.4 Sampling 

We would like to estimate the statistical error of Z*. Ideally, we want the standard deviation 
of the estimate smaller than some given tolerance. The major contribution to the sampling er- 
ror of Z* comes from that of (VF)j, because of the cancelation and subsequent multiplication 
by e~ 2 . Thus, we neglect the statistical error of Dj in estimating the error of gj. So then, the 
statistical error of Z* comes from the sample average of Aj = PDjfflF)™, n = 1, 2, . . . , N, 
where N is the sample size. The statistical error is defined by (maxo<j<j error bar of Aj)r 2 , 
where an error bar is an estimate of 1 standard deviation. Such an estimate can be obtained 
using block averaging as in Ref. [10] (Appendix D.3). In general, 32 blocks is a reasonable 
choice. 

At each iteration, the configuration x from the previous iteration could be used to start 
the equilibration of the molecular dynamics. Thus, it is necessary that values of x be stored 
such that = Zj, j = 0, 1, . . . , J. It is reasonable to expect less equilibration time is 
needed in later iterations as the path converges. 

5 Numerical tests 
5.1 An artificial problem 

As an example to illustrate our method, consider a problem finding the MFTP and MFEP 
for the potential energy function 

U(x,y) = -4exp(-4x 2 -(y-2.75) 2 )-5exp(-(a;-l) 2 -(y-0.15) 2 ) 

- 5exp(-(x + l) 2 -y 2 ) + 8exp(-x 2 - (y + 0.5) 2 ) + 0.001(x 4 + y 4 ) 

where the energy unit is kcal/mol and the mass matrix M has identical diagonal entries. Un- 
less specifically mentioned, the inverse temperature = 0.59595 kcal/mol, corresponding 
to 300 K. In particularly, we take collective variables ( = = (x,y). In this case, the 

MFEP becomes a minimum energy path (MEP). Alternatively, an MEP can be considered 
as an MFEP for which we have an accurate estimate of -F(C)- 

In [3j we show an MEP connecting two local minima through the third local mini- 
mum. The MEP has a cusp at the intermediate minimum. The MEPs are computed 
using the simplified string method with piecewise linear interpolation and equal arc length 
normalization. The time step r 2 = 0.01. The iteration is stopped if d < 0.005, where 
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d = max <j<j t~ 2 \ZJ cw — Zj\. From the figure, we can see that the cusp is missing if the 
number of images (J + 1 = 10) is too small. Also, the MEP does not go through the 
intermediate local minimum as it should, even with many images ( J + 1 = 80). 




Figure 3: Minimum energy path obtained using the simplified string method. The initial 
path is the straight line between (—1,0) and (1,0). The path is discretized into J+l images. 
Four figures are generated using J + 1 = 10, 20, 40, 80 images, respectively. 



A calculation (not shown here) similar to that for [3] was done for the MFTP. The MFTP 
is calculated using the semi-implicit simplified string method described in Sec. |4j The MFTP 
can be resolved using a relatively small set of images, for example, the MFTP calculated by 
only 10 images (J = 9) is almost indistinguishable from the one calculated using 80 images 
(J = 79). The MFTP avoids the cusp problem. 

The MFTP generates different paths at different temperature. [4] shows MFTPs at 3 K, 
30 K, 300 K, 3000 K, 30000 K, respectively. It is clear that the MFTP is close to the MEP at 
low temperature (3K) and is close to a straight line at high temperature (30000 K), which 
is what we expect. 

An FTS method path is expected to be similar to an MFEP for this example. 
5.2 Phi, psi for alanine dipeptide in vacuum 

For comparison with the MFEP, we study alanine dipeptide at 300 K in vacuum [20]. We 
compare the MFEP and the MFTP with two dihedral angles <fi and if) as collective variables. 
All simulations were performed using the CHARMM simulation program [JJ [3] and the full- 
atom representation of the molecule in the CHARMM force field [TSJ HHI - Langevin dynamics 
with friction coefficient 10.0 ps _1 and time step 1.0 fs was used. For the calculation of V-F 



16 



Figure 4: Maximum flux transition path obtained using the semi-implicit simplified string 
method. Here we used the same initial path and the same stopping criterion for convergence 
as for|3} The MFTPs are generated using 20 images at 3 K, 30 K, 300 K, 3000 K, 30000 K, re- 
spectively (which roughly correspond to = 0.006,0.06,0.6,6,60 kcal/mol.) The contour 
lines are separated by 0.25 kcal/mol. 



and D, harmonic potentials as in Eq. (14) were added involving the dihedral angles <fi and 
i/j with force constant 1000 kcal/(mol rad 2 ) (corresponding to e = 1°). 

The initial path in collective variable space is a straight line between two points in </> — ip 
space. The path is discretized into J + 1 images. The configuration of alanine dipeptide at 
each image along the initial path is built using the IC module in CHARMM with dihedral 
angles fixed at the interpolated values. Then follow 1000 steps of minimization and 50,000 
steps of heating before the iteration starts. Each iteration of the path involves 50,000 steps of 
equilibration and 500,000 steps of sampling. The configuration at the final step of sampling 
in the previous iteration is used as the initial configuration for the equilibration in the next 
iteration. 

We begin by comparing the MFTP and MFEP from C 7eq to C 7ax . The MFEP is calculated 
using the simplified string method with linear interpolation between images and equal arc 
length normalization. The MFTP is calculated using the semi-implicit simplified string 
method. In [5j the initial path is the straight line between (—83.2°, 74.5°) and (70°,— 70°), 
which were determined as Cj cq and CV ax in Ref . (20] . The path is discretized into 20 images. 
The time step r 2 = 0.16 in CHARMM time units squared, or r 2 = (19.56 fs) 2 . The statistical 
error estimated by block averaging using 32 blocks is ±0.00577°. The iteration is stopped 
if d < 0.02. (The tolerance value should be chosen properly since the statistical error will 
eventually dominate the other errors so that d fluctuates about a positive number.) It takes 
34 and 31 iterations to converge for the MFTP and MFEP, respectively. The computational 
cost for two methods is comparable. The path calculated for this problem by the FTS 
method using the CHARMM force field is given in Figure 5 of Ref. [29J . 
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Figure 5: Maximum flux transition path and minimum free energy path from CV cq to CV ax for 
alanine dipeptide in vacuum at 300 K. Triangles are images of the initial path; rectangles are 
the images of the maximum flux transition path; and circles are the images of the minimum 
free energy path. The contours are those for the zero-temperature free energy (adiabatic 
energy). The contour lines are separated by 0.6kcal/mol. 



Next we compare the MFTP and MFEP from CV cq to C' 7 . In particular, we calculate 
the transition path Cj eq -C 7ax -C' 7eq , in which C 7ax serves as an intermediate metastable state. 
The initial path is taken to be the straight line between (-80°, 80°) and (190°, -190°). g 
shows the MFTP and MFEP generated using 40 images. The time step r 2 = (19.56 fs) 2 . 
The iteration is stopped if d < 0.02. It takes 35 and 44 iterations for the MFTP and MFEP 
to converge, respectively. It is evident that the MFTP is more efficient than the MFEP in 
this case. 

5.3 Phi, psi for alanine dipeptide in solution 

We also test our method for alanine dipeptide solvated in explicit water. Again, the backbone 
dihedrals <fi and ip are used as collective variables to describe the transition. The initial 
paths are straight lines connecting two points among (—77°, 138°), (55°, 48°), (60°,— 72°), 
and (-77°, -39°) in 0,^)-space. 

For preparing the simulation, each starting structure for alanine dipeptide with con- 
strained and ip angels is solvated in a (20 x 18 x 15) A 3 box with 191 TIP3 [15] water 
molecules and equilibrated for 50,000 ps. The molecular dynamics are carried out with the 
CHARMM program under the CHARMM22 force field. Periodic boundary conditions are 
used and the electrostatic interactions are treated with the particle mesh Ewald method [U] . 
The system is simulated at a constant pressure of 1.0 atm and a constant temperature 300 K 
with the algorithm based on Hoover's methods. We use a 1-fs timestep with the SHAKE [30] 
algorithm to keep all bonds involving hydrogen atoms at fixed lengths. 



18 




-100 -50 50 100 150 200 



4> 

Figure 6: Maximum flux transition path and minimum free energy path for alanine dipeptide 
from CV eq to C' 7eq passing by CV ax in vacuum at 300 K. The figure is generated using 40 images. 
Triangles are the images for the initial path; rectangles are the images of the maximum flux 
transition path; and circles are the images of the minimum free energy path. The contours are 
those for the zero-temperature free energy. The contour lines are separated by 0.6kcal/mol. 



In [7j four MFTPs are calculated using the semi-implicit simplified string method. Each 
iteration involves 50,000 steps of equilibration and 500,000 steps of sampling. The transition 
paths are the result of 50 cycles of iteration. The path C-j eq -oiL calculated for this problem 
by the FTP method using the CHARMM force field is given in Figure 12 of Ref . [29] . The 
MFTP is similar to the FTS method path. 
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Figure 7: Maximum flux transition paths for alanine dipeptide in solution. The transition 
paths are calculated by the semi-implicit simplified string method with the nearby straight 
lines as initial paths. 

A Derivation of Brownian dynamics approximation 

The quasi-committor is related to a full phase-space committor in Ref. [33] (Sec. 6. 2) q* 
defined as follows: 

q*(x,p) = Pr(X(t) reaches B before A | X(0) = x, P(0) = p). 

Note that q*(x,p) = or 1, because the dynamical equation is deterministic. By definition, 
the quasi-committor q(Q = (q*(x,p))^. 

It is not difficult to show that q{£{x)) approximates q*{x,p) in the sense that it mini- 
mizes (\q(£(x)) — q*(x,p)\ 2 ) over all q(()- However, this is not useful for determining q(() 
because q*(x,p) is too costly to compute. On the other hand, it is possible to find a best 
approximation to q*(x,p) in another sense. Because q* is constant on a trajectory, we have 

= ^q*(X(t),P(t)) = (Lq*)(X(t),P(t)) where L = {M~ l p) -V x -U x - V p . 

Consequently, q* satisfies the stationary Liouville equation 

Lq* =0, q* = on A, q* = 1 on B. 

Since we do know Lq* = 0, we seek instead an approximation q that minimizes I(q) = 
(\L(q(£(x)) — q*{x,p))\ 2 ) , a standard tactic in numerical analysis. As shown in Sec. III.B of 
Ref. [20], this simplifies to 

/(g) = ^(|M- 1 / 2 V a; g(^))| 2 ), 
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which is to be as small as possible. A low value for I(q) is attained by having q(() increase 
monotonically from the value on to the value 1 on B^, which is consistent with the 
prescription given earlier that £(x) be chosen so that the quasi-committor has no local 
minima or maxima outside of and B^. 

The functional I(q) can be expressed in terms of collective variables ( as given by Eq. ^ 
and shown in in Eq. (15) of Ref. [2D] . 



B Derivation of lasting hitting-point distribution 

The proof of Proposition 5 in Ref. [H] (p. 158) analyzes the flux of reactive trajectories. The 
flux J(Q gives the rate at which such trajectories cross an arbitrary surface E that separates 
collective variable space into two parts, one containing A^ and the other containing B^ via 
the integral J 2 J(() -h(()dS^ where h(Q points to the side containing B^. The proof actually 
looks not at all crossings but only those occurring within a vanishingly small time interval 
before the last crossing — see Eq. (50) of Ref. [8]. Therefore, it considers the net flux only in 
this limiting sense. As the length of the time interval r — > 0, the positions of these crossings 
all converge to the position of the last crossing. So, indeed, one gets the flux of the last 
hitting point from Proposition 5 of Ref. [8]. The result given in Ref. [E] (Eq. (39)), as well as 
in Ref. [2T] (Eq. (6), Eq. (A12)), and Ref. [33] (Eq. (62)), is that the last hitting-point flux 
for reactive trajectories is J(() = p^(C)-D(C)Vg(C)- (Proposition 4 of Ref. [8] does not apply 
to the infinitely damped case of Langevin dynamics.) That the expression for J(Q also gives 
the net flux of reactive trajectories is Eq. (32) of Ref. [33]. Also, the formula for in 



Sec. 3^2 agrees in the special case D = I with that for the first hitting point distribution 
given in Ref. [37] (Appendix B). Last and first are the same for reversible dynamics like 
Brownian dynamics. Finally, there is an example in Metzner, Schutte, and Vanden-Eijnden 
(2006) section HI.C, where it is suggested to use h ■ J. 



C Derivation of the maximum flux condition 

We have from (jl]) and ^ that the normal flux is 

j(C) = CO n^e-^^g s ( ( r(C))V ( 7(C) T ^(C)V ( 7(C)/|V ( 7(C)| (15) 

where a(() is defined implicitly by q(() = q(a(Q). And for each point Z(s) on the desired 
path, the condition to be satisfied ^ is Vj(Z(s)) || n(s). We also approximate -D(C) by 

D(a(Q) where D{s) == D(Z(s)). Furthermore, the assumption (8) that isocommittors are 
planar implies 

«))-(C-«))) = 0. (16) 



Differentiating Eq. (16) w.r.t. we get 



(n» • (C - Z{a)) - n{a) ■ Z s {a))Va + n{a) = 0, 
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where the argument ( of a has been omitted, whence 

Va = (n(a) ■ Z s (a) - n s (a) ■ (C - Z{a)))- 1 n{(j). 



(17) 



Substituting Eq. (17) into Eq. (15) and replacing D(() by D(a(()) } the normal flux becomes 
where 

tp(s, C) = const^ F ^q s (s)(n(s) ■ Z s {s) - n s (s) ■ (( - Z(s)y l n(s) T D(s)n(s)/\n(s)\. 
Note that 



V c <^ 



-(3VF 



and 



Thus, we have 



V c <p 



n ■ Z s - n s ■ (C - Z) ' 



-PVF{Z) + 



n T Z x 



and 



3 



V9 s (s,Z)n 



-/3VF(Z) + + 



n T Z s ip(s,Z)n(s) T Z s (s)' 



Hence, the condition is that 
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Captions 

Figure 1: A schematic illustration of a poor choice of collective variables. The horizontal 
axis is collective variables, and the vertical axis is unrepresented degrees of freedom. The 
collective variables fail to indicate the progress of the reaction. 

Figure 2: Shading indicates contours of free energy, thin curves denote isocommittors, 
ellipses enclose concentrations of crossing points from reactive trajectories, and the thick 
curve is the center. 
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Figure 3: Minimum energy path obtained using the simplified string method. The initial 
path is the straight line between (—1, 0) and (1,0). The path is discretized into J+l images. 
Four figures are generated using J + 1 = 10, 20, 40, 80 images, respectively. 

Figure 4: Maximum flux transition path obtained using the semi-implicit simplified string 
method. Here we used the same initial path and the same stopping criterion for convergence 
as for § The MFTPs are generated using 20 images at 3 K, 30 K, 300 K, 3000 K, 30000 K, re- 
spectively (which roughly correspond to f3~ l = 0.006,0.06,0.6,6,60 kcal/mol.) The contour 
lines are separated by 0.25 kcal/mol. 

Figure 5: Maximum flux transition path and minimum free energy path from CV cq to 
C*7 ax for alanine dipeptide in vacuum at 300 K. Triangles are images of the initial path; 
rectangles are the images of the maximum flux transition path; and circles are the images of 
the minimum free energy path. The contours are those for the zero-temperature free energy 
(adiabatic energy). The contour lines are separated by 0.6 kcal/mol. 

Figure 6: Maximum flux transition path and minimum free energy path for alanine 
dipeptide from Cj cq to C' 7cq passing by CV ax in vacuum at 300 K. The figure is generated 
using 40 images. Triangles are the images for the initial path; rectangles are the images of 
the maximum flux transition path; and circles are the images of the minimum free energy 
path. The contours are those for the zero-temperature free energy. The contour lines are 
separated by 0.6 kcal/mol. 

Figure 7: Maximum flux transition paths for alanine dipeptide in solution. The transition 
paths are calculated by the semi-implicit simplified string method with the nearby straight 
lines as initial paths. 
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